#################################################################
############Analysis file for Objectified & Dehumanized###########
##########################Nov. 21, 2020###########################
##################################################################

library(stargazer)
library(effects)
library(jtools)
library(ATE)
library(readstata13)
library(table1)
library(ggplot2)
library(interplot)
library(sjPlot)
library(sjmisc)
library(ggpubr)
library(interactions)

#Import Cleaned Data

ObjectData <- read.dta13("O&D, clean.dta")
head(ObjectData)


#Demographic Table 

demo <- subset(ObjectData, select=c("age", "conservative", "genderfemale", "education", "income", "white", "latinx", "evangelical", "Democrat", "Republican"))

summary(demo)

demo$genderfemale <-
  factor(demo$genderfemale, levels=c(1,0), 
         labels=c("Women",
                  "Men"))

demo$white <- 
  factor(demo$white, levels=c(1,0), 
         labels=c("White", 
                  "Non-White"))

demo$latinx <- 
  factor(demo$latinx, levels=c(1,0), 
         labels=c("Latinx", 
                  "Non-Latinx"))

demo$education <- 
  factor(demo$education, levels=c(1,2,3,4,5,6), 
         labels=c("No H.S. diploma", 
                  "H.S. diploma",
                  "Some college",
                  "Associate's degree",
                  "Bachelor's degree",
                  "Graduate degree"))

demo$income <- 
  factor(demo$income, levels=c(1,2,3,4,5,6,7), 
         labels=c("Below $20,000", 
                  "$20,001 to $40,000",
                  "$40,001 to $60,000",
                  "$60,001 to $80,000",
                  "$80,001 to $100,000",
                  "$100,001 to $120,00",
                  "Over $120,000"))

demo$evangelical <- 
  factor(demo$evangelical, levels=c(1,0), 
         labels=c("Evangelical", 
                  "Non-Evangelical"))

demo$Democrat <- 
  factor(demo$Democrat, levels=c(1,0), 
         labels=c("Democrat", 
                  "Republican or Independent"))

demo$Republican <- 
  factor(demo$Republican, levels=c(1,0), 
         labels=c("Republican", 
                  "Democrat or Independent"))

label(demo$genderfemale) <- "Gender"
label(demo$white) <- "Race"
label(demo$latinx) <- "Latinx"
label(demo$education) <- "Education"
label(demo$income) <- "Income"
label(demo$evangelical) <- "Evangelical Identification"
label(demo$Democrat) <- "Democrat"
label(demo$Republican) <- "Republican"


table <- table1(~ genderfemale + white + latinx + education + income + evangelical + Democrat + Republican + age, data=demo, output="html", export="table1.html")

table

###Set treatment and genderfemale as factors 

is.numeric(ObjectData$treatment)
is.numeric(ObjectData$genderfemale)

ObjectData$treatment.f <- factor(ObjectData$treatment, labels = c("control", "treatment"))
head(ObjectData)
summary(ObjectData$treatment.f)

ObjectData$genderfemale.f <- factor(ObjectData$genderfemale, labels = c("Men", "Women"))
head(ObjectData)
summary(ObjectData$genderfemale.f)

#Difference in support for women in politics (summated rating scale)

M1_H1 <- lm(WIP ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
M2_H1 <- lm(WIP ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f, data=ObjectData)
  
summary(M1_H1)
summary(M2_H1)  


#Regression Table #1 

stargazer(M1_H1, M2_H1, type="text", align=TRUE, title="Results",
          dep.var.labels=c("WIP Support"),
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "TreatmentXGender"),
          out="Hyp1_Table1.htm")

M1_H1_evan <- lm(WIP ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + evangelical, data=ObjectData)
M2_H1_evan <- lm(WIP ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f + evangelical, data=ObjectData)

stargazer(M1_H1_evan, M2_H1_evan, type="text", align=TRUE, title="Results",
          dep.var.labels=c("WIP Support"),
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "Evangelical", "TreatmentXGender"),
          out="Hyp1_Table1_evangelical.htm")

#Visualization #1 and #2

Viz1 <- effect_plot(M1_H1, pred = "treatment.f", interval = TRUE, y.label = "Support for WIP (1-5)", x.label = "Treatment", pred.labels = c("Objectifying Images", "Neutral Images"), force.cat = TRUE, point.shape = TRUE)+
        ylim(3,4)
ggsave("Plot1.tiff", width = 4.5, height = 3.5, dpi=700)


Viz2 <- cat_plot(M2_H1, pred = treatment.f, modx = genderfemale.f, y.label = "Support for WIP (1-5)", x.label = "Treatment", legend.main="Gender")

Viz2

ggsave("Plot2.tiff", width = 4.5, height = 3.5, dpi=700)

#Difference in evaluation of women candidates (summated rating scale)

M1_H2 <- lm(WPE ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
M2_H2 <- lm(WPE ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f, data=ObjectData)

summary(M1_H2)
summary(M2_H2)  

#Regression Table #2

stargazer(M1_H2, M2_H2, type="text", align=TRUE, title="Results",
          dep.var.labels=c("Women Politicians Evaluation"),
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "TreatmentXGender"),
          out="Hyp2_Table1.htm")

M1_H2_evan <- lm(WPE ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + evangelical, data=ObjectData)
M2_H2_evan <- lm(WPE ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f + evangelical, data=ObjectData)

stargazer(M1_H2_evan, M2_H2_evan, type="text", align=TRUE, title="Results",
          dep.var.labels=c("Women Politician Evaluation"),
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "Evangelical", "TreatmentXGender"),
          out="Hyp2_Table2_evangelical.htm")

#Visualization #3 and #4

Viz3 <- effect_plot(M1_H2, pred = "treatment.f", interval = TRUE, y.label = "Evaluations of Women Politicians (0-10)", x.label = "Treatment", force.cat = TRUE, point.shape = TRUE) + 
  ylim(3.5,5)
Viz3
ggsave("Plot3.tiff", width = 4.5, height = 3.5, dpi=700)


Viz4 <- cat_plot(M2_H2, pred = treatment.f, modx = genderfemale.f, y.label = "Evaluations of Women Politicians (0-10)", x.label = "Treatment", legend.main="Gender")
ggsave("Plot4.tiff", width = 4.5, height = 3.5, dpi=700)


#Dehumanization of women candidates (summated rating scale)

M1_H3 <- lm(dehumanization ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
M2_H3 <- lm(dehumanization ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f, data=ObjectData)

summary(M1_H3)
summary(M1_H3)


#Regression Table #3

stargazer(M1_H3, M2_H3, type="text", align=TRUE, title="Results",
          dep.var.labels=c("Overall Dehumanization"),
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "TreatmentXGender"),
          out="Hyp3_Table1.htm")

M1_H3_evan <- lm(dehumanization ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + evangelical, data=ObjectData)
M2_H3_evan <- lm(dehumanization ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f + evangelical, data=ObjectData)

stargazer(M1_H3_evan, M2_H3_evan, type="text", align=TRUE, title="Results",
          dep.var.labels=c("Overall Dehumanization"),
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "Evangelical", "TreatmentXGender"),
          out="Hyp3_Table2_evangelical.htm")

#Visualization #5 and #6

Viz5 <- effect_plot(M1_H3, pred = "treatment.f", interval = TRUE, y.label = "Overall Dehumanization of Women (1-4)", x.label = "Treatment", force.cat = TRUE, point.shape = TRUE)+
  ylim(2.25,2.75)
Viz4
ggsave("Plot5.tiff", width = 4.5, height = 3.5, dpi=700)


Viz6 <- cat_plot(M2_H3, pred = treatment.f, modx = genderfemale.f, y.label = "Overall Dehumanization of Women (1-4)", x.label = "Treatment", legend.main="Gender")
ggsave("Plot6.tiff", width = 4.5, height = 3.5, dpi=700)


########################################################            
#####Ordered Logistic Regression- Robustness Checks#####
########################################################

##all variables are continuous 

#Additional Analyses- Individual candidates 

Warren <- lm(Warren ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
Harris <- lm(Harris ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
Haley <- lm(Haley ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
AOC <- lm(AOC ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
Pelosi <- lm(Pelosi ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)


stargazer(Warren, Harris, Haley, AOC, Pelosi, type="text", align=TRUE, title="Results",
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption"),
          out="Indiv_Cans.htm")

#Dehumanization- Individual Scales 

Animalistic_1 <- lm(animalistic ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
Animalistic_2 <- lm(animalistic ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f, data=ObjectData)

Mechanistic_1 <- lm(mechanistic ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume, data=ObjectData)
Mechanistic_2 <- lm(mechanistic ~ treatment.f + genderfemale.f + age + white + education + Republican + mediaconsume + treatment.f*genderfemale.f, data=ObjectData)

stargazer(Animalistic_1, Animalistic_2, Mechanistic_1, Mechanistic_2, type="text", align=TRUE, title="Results",
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "TreatmentXGender"),
          out="Indiv_Dehumanization.htm")


######Combine Regression Tables#######

stargazer(M1_H1, M2_H1, M1_H2, M2_H2, M1_H3, M2_H3, type="text", align=TRUE, title="Results",
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "TreatmentXGender"),
          out="Combined_Main_Table.htm")

####Combine Regression Tables with Evangelical Identification######


stargazer(M1_H1_evan, M2_H1_evan, M1_H2_evan, M2_H2_evan, M1_H3_evan, M2_H3_evan, type="text", align=TRUE, title="Results",
          covariate.labels=c("Treatment", "Gender", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "Evangelical", "TreatmentXGender"),
          out="Combined__Table_Evangelical.htm")

###################Combined Figures########################


Combined_Fig1 <- ggarrange(Viz1, Viz3, Viz5, 
          common.legend = TRUE, legend = "bottom")

ggsave("Combined_Fig1.tiff", width = 7, height = 8, dpi=700)


Combined_Fig2 <- ggarrange(Viz2, Viz4, Viz6, 
                           common.legend = TRUE, legend = "bottom")

ggsave("Combined_Fig2.tiff", width = 7, height = 8, dpi=700)

######Continuous Measures of Gender Identity 

M1_H1_cont <- lm(WIP ~ treatment.f + genderfemale.f + femininity + masculinity + age + white + education + Republican + mediaconsume, data=ObjectData)
M2_H1_cont <- lm(WIP ~ treatment.f + genderfemale.f + femininity + masculinity + age + white + education + Republican + mediaconsume + treatment.f*femininity + treatment.f*masculinity, data=ObjectData)

M1_H2_cont <- lm(WPE ~ treatment.f + genderfemale.f + femininity + masculinity + age + white + education + Republican + mediaconsume, data=ObjectData)
M2_H2_cont <- lm(WPE ~ treatment.f + genderfemale.f + femininity + masculinity + age + white + education + Republican + mediaconsume + treatment.f*femininity + treatment.f*masculinity, data=ObjectData)

M1_H3_cont <- lm(dehumanization ~ treatment.f + genderfemale.f + femininity + masculinity + age + white + education + Republican + mediaconsume, data=ObjectData)
M2_H3_cont <- lm(dehumanization ~ treatment.f + genderfemale.f + femininity + masculinity + age + white + education + Republican + mediaconsume + treatment.f*femininity + treatment.f*masculinity, data=ObjectData)


stargazer(M1_H1_cont, M2_H1_cont, M1_H2_cont, M2_H2_cont, M1_H3_cont, M2_H3_cont, type="text", align=TRUE, title="Results",
          covariate.labels=c("Treatment", "Gender", "Fem Scale", "Masc Scale", "Age", "Race", "Education", "Republican", 
                             "Media Consumption", "TreatmentXFem", "TreatmentXMasc"),
          out="Combined_Main_Table_Cont.htm")



####bivariate correlations 

#by gender and treatment group, between the dehumanization scales, WIP, evaluations 

cor(ObjectData$WIP, ObjectData$dehumanization, method = c("pearson"))
cor(ObjectData$WIP, ObjectData$dehumanization,use = "complete.obs", method = c("pearson"))
cor(ObjectData$genderfemale, ObjectData$treatment,use = "complete.obs", method = c("pearson"))
cor(ObjectData$genderfemale, ObjectData$treatment, use = "complete.obs", method = c("pearson"))



